# Codes for store-group price index aggregation to store-level 
# JHL 

# Time printing function 
verbose.Print.Time <- function(verbose, msg=NA, prev.time=NA) {
    cur.time <- proc.time()[3]
    if (verbose) {
        if (is.na(prev.time)) { # no time to print w/o prev.time 
            cat(sprintf("%s", msg))
        } else {
            if (is.na(msg)) { # if no message, end of timing, just print time
                cat(sprintf("(%1.0fs), ", cur.time-prev.time))
            } else { # usual case, print time since last message and message
                cat(sprintf("(%1.0fs); %s ", cur.time-prev.time, msg))
            }
        }	
    }
    return(cur.time)
}

# Input store x product group indices, calculate sales by store, select only stores that are balanced by group and have observations every period 
PI.pre <- function(PI) {

  # Eliminate extra variables if necessary
  PI[, x:= NULL]
  PI[, group:=NULL]
   
	PI[, VUL_t_1 := c(0,VUL_t[-.N]), by = Group] 
	# Sum first within each unit (store) because next step removes groups that are unbalanced
	PI[,`:=`(Sum_TVC_t = sum(TVC_t), Sum_TVC_t_1 = sum(TVC_t_1),
	         Sum_VUL_t = sum(VUL_t), Sum_VUL_t_1 = sum(VUL_t_1),
			 Sum_P_y_lm = sum(P_y_lm), Sum_P_y_glm = sum(P_y_glm),
			 Sum_PI_start = sum(PI_start), Sum_PI_end = sum(PI_end),
			 Sum_PI_gl_start = sum(PI_gl_start), Sum_PI_gl_end = sum(PI_gl_end),
			 Sum_P_y_fix_lm = sum(P_y_fix_lm)), by = quarter] 

	PI_no_nan <- PI
	
	return(PI_no_nan) 
}

# Calculate VUL to TVC ratio if called, by store 
Calc.TVC.VUL <- function(PI) {
	if ( !is.null(PI) ) {
		TVC_VUL <- PI[, list(quarter, Sum_TVC_t, Sum_VUL_t)][1:40]
		PI[, Sum_VUL_t := NULL]
		setnames(TVC_VUL, c("Sum_TVC_t", "Sum_VUL_t"), 
						  c("TVC_t", "VUL_t"))
		 
		Total <- TVC_VUL[, list(TVC_t = sum(TVC_t), 
								VUL_t = sum(VUL_t))]
		Total[, quarter := 0L]
		setcolorder(Total, c("quarter", "TVC_t","VUL_t"))
		l = list(Total, TVC_VUL)
		rm(Total)

		TVC_VUL <- rbindlist(l, use.names = TRUE)
		rm(l)
		TVC_VUL[, `:=`(Ratio_TVC = VUL_t/TVC_t)]
		return(TVC_VUL)
	}
}

# Construct multiple indices: yes index is Tornqvist, nsub is Laspeyres 	
# Construct yes indices 
Calc.Tornqvist <- function(PI_no_nan) {
    if ( !is.null(PI_no_nan) ) {       

		# 1. Laspeyres construction 
		PI_l_no_nan <- PI_no_nan[P_t != "NA",] # Keep only if index was calculated 
		# Keep only Groups that are appear in entire sample
		PI_l_no_nan[,max := max(quarter), by = Group]
		PI_l_no_nan <- PI_l_no_nan[max == 40, ]
		PI_l_no_nan[, max:=NULL]	

		if(nrow(PI_l_no_nan)>0){
			PI_l_no_nan[, `:=` (S_TVC_t = (TVC_t/Sum_TVC_t), S_TVC_t_1 = (TVC_t_1/Sum_TVC_t_1),
					   S_VUL_t = (VUL_t/Sum_VUL_t), S_VUL_t_1 = (VUL_t_1/Sum_VUL_t_1)),
					   by = quarter]	
			
			PI_l_no_nan[, `:=` (growth = 1.0, exp_TVC = 0.0,
							  exp_VUL = 0.0), by = quarter]

			PI_l_no_nan[quarter >1, `:=` (growth = P_t/P_t_1,
					  exp_TVC = (S_TVC_t+S_TVC_t_1)/2, exp_VUL = (S_VUL_t+S_VUL_t_1)/2),
					   by = quarter]
			
			# Each term with exponent using average shares 
			PI_l_no_nan[, `:=`(pTornTVC_t = growth^exp_TVC,
							 pTornVUL_t = growth^exp_VUL),
							 by = quarter]
			
			# Log, sum, then unlog (instead of using product)
			PI_l_no_nan[, `:=`(logpT_TVC = log(pTornTVC_t),
							 logpT_VUL = log(pTornVUL_t) )]
			
			# Split out and merge back in after ratios are calculated
			PI_agg <- PI_l_no_nan[, list(PT_TVC = sum(logpT_TVC),
									   PT_VUL = sum(logpT_VUL)),
									   by = quarter]

			PI_l_no_nan <- PI_l_no_nan[, list(quarter, Group, TVC_t, VUL_t, Sum_TVC_t, Sum_TVC_t_1, 
										  Sum_VUL_t, Sum_VUL_t_1, Sum_P_y_lm, Sum_P_y_glm, Sum_PI_start, Sum_PI_end, Sum_PI_gl_start, Sum_PI_gl_end, Sum_P_y_fix_lm )]
				  
			# Calculate ratio of revenue included into the index construction after dropping balanced groups / selection criteria in PI.pre
			PI_l_no_nan[, ':=' (Ratio_TVC = sum(VUL_t)/Sum_TVC_t, Ratio_VUL = sum(VUL_t)/Sum_VUL_t), by = quarter]

			PI_agg[, `:=`(PT_TVC = exp(PT_TVC), PT_VUL = exp(PT_VUL) )] # logged, can use sum, then unlog, for products 
			
			# Create yes indices, normalize q1 to 1, then multiply each PT_yes_TVC by the last term in the equation
			for (i in 1:nrow(PI_agg)) {
				PI_agg[i, `:=`(PT_yes_TVC = ifelse(quarter == 1, 1.0, 
							   PI_agg[i, PT_TVC]*PI_agg[i-1, PT_yes_TVC]),
							   PT_yes_VUL = ifelse(quarter == 1, 1.0, 
							   PI_agg[i, PT_VUL]*PI_agg[i-1, PT_yes_VUL]) )]

			}

			setkey(PI_agg,quarter)
			setkey(PI_l_no_nan, quarter)
		}
		
		# 2. Geometric Laspeyres construction 
		PI_gl_no_nan <- PI_no_nan[P_gl_t != "NA",] # Keep only if index was calculated 
		# Keep only Groups that are appear in entire sample
		PI_gl_no_nan[,max := max(quarter), by = Group]
		PI_gl_no_nan <- PI_gl_no_nan[max == 40, ]
		PI_gl_no_nan[, max:=NULL]		
		
		if (nrow(PI_gl_no_nan)>0){
			PI_gl_no_nan[, `:=` (S_TVC_t = (TVC_t/Sum_TVC_t), S_TVC_t_1 = (TVC_t_1/Sum_TVC_t_1),
					   S_VUL_t = (VUL_t/Sum_VUL_t), S_VUL_t_1 = (VUL_t_1/Sum_VUL_t_1)),
					   by = quarter]	
			
			PI_gl_no_nan[, `:=` (growth = 1.0, exp_TVC = 0.0,
							  exp_VUL = 0.0), by = quarter]

			PI_gl_no_nan[quarter >1, `:=` (growth_gl = P_gl_t/P_gl_t_1,
					  exp_TVC = (S_TVC_t+S_TVC_t_1)/2, exp_VUL = (S_VUL_t+S_VUL_t_1)/2),
					   by = quarter]
			
			# Each term with exponent using average shares 
			PI_gl_no_nan[, `:=`(p_gl_TornTVC_t = growth_gl^exp_TVC,
							 p_gl_TornVUL_t = growth_gl^exp_VUL),
							 by = quarter]
			
			# Log, sum, then unlog (instead of using product)
			PI_gl_no_nan[, `:=`(logpglT_TVC = log(p_gl_TornTVC_t),
							 logpglT_VUL = log(p_gl_TornVUL_t) )]
			
			# Split out and merge back in after ratios are calculated
			PI_gl_agg <- PI_gl_no_nan[, list(PglT_TVC = sum(logpglT_TVC),
									   PglT_VUL = sum(logpglT_VUL)),
									   by = quarter]

			PI_gl_no_nan <- PI_gl_no_nan[, list(quarter, Group, TVC_t, VUL_t, Sum_TVC_t, Sum_TVC_t_1, 
										  Sum_VUL_t, Sum_VUL_t_1, Sum_P_y_lm, Sum_P_y_glm, Sum_PI_start, Sum_PI_end, Sum_PI_gl_start, Sum_PI_gl_end, Sum_P_y_fix_lm )]
				  
			# Calculate ratio of revenue included into the index construction after dropping balanced groups / selection criteria in PI.pre
			PI_gl_no_nan[, ':=' (Ratio_TVC = sum(VUL_t)/Sum_TVC_t, Ratio_VUL = sum(VUL_t)/Sum_VUL_t), by = quarter]

			PI_gl_agg[, `:=`(PglT_TVC = exp(PglT_TVC), PglT_VUL = exp(PglT_VUL) )] # logged, can use sum, then unlog, for products 
			
			# Create yes indices, normalize q1 to 1, then multiply each PT_yes_TVC by the last term in the equation
			for (i in 1:nrow(PI_gl_agg)) {
				PI_gl_agg[i, `:=`(PglT_yes_TVC = ifelse(quarter == 1, 1.0, 
							   PI_gl_agg[i, PglT_TVC]*PI_gl_agg[i-1, PglT_yes_TVC]),
							   PglT_yes_VUL = ifelse(quarter == 1, 1.0, 
							   PI_gl_agg[i, PglT_VUL]*PI_gl_agg[i-1, PglT_yes_VUL]) )]

			}

			setkey(PI_gl_agg,quarter)
			setkey(PI_gl_no_nan, quarter)
		}
		
		# 3. Fixed index construction 
		PI_fix_no_nan <- PI_no_nan[P_fix_t != "NA",] # Keep only if index was calculated 
		# Keep only Groups that are appear in entire sample
		PI_fix_no_nan[,max := max(quarter), by = Group]
		PI_fix_no_nan <- PI_fix_no_nan[max == 40, ]
		PI_fix_no_nan[, max:=NULL]		
		
		if (nrow(PI_fix_no_nan)!=0) {
			PI_fix_no_nan[, `:=` (S_TVC_t = (TVC_t/Sum_TVC_t), S_TVC_t_1 = (TVC_t_1/Sum_TVC_t_1),
					   S_VUL_t = (VUL_t/Sum_VUL_t), S_VUL_t_1 = (VUL_t_1/Sum_VUL_t_1)),
					   by = quarter]	
			
			PI_fix_no_nan[, `:=` (growth = 1.0, exp_TVC = 0.0,
							  exp_VUL = 0.0), by = quarter]

			PI_fix_no_nan[quarter >1, `:=` (growth_fix = P_fix_t/P_fix_t_1,
					  exp_TVC = (S_TVC_t+S_TVC_t_1)/2, exp_VUL = (S_VUL_t+S_VUL_t_1)/2),
					   by = quarter]
		
			# Each term with exponent using average shares 
			PI_fix_no_nan[, `:=`(p_fix_TornTVC_t = growth_fix^exp_TVC,
							 p_fix_TornVUL_t = growth_fix^exp_VUL),
							 by = quarter]
			
			# Log, sum, then unlog (instead of using product)
			PI_fix_no_nan[, `:=`(logpfixT_TVC = log(p_fix_TornTVC_t),
							 logpfixT_VUL = log(p_fix_TornVUL_t) )]
			
			# Split out and merge back in after ratios are calculated
			PI_fix_agg <- PI_fix_no_nan[, list(PfixT_TVC = sum(logpfixT_TVC),
									   PfixT_VUL = sum(logpfixT_VUL)),
									   by = quarter]

			PI_fix_no_nan <- PI_fix_no_nan[, list(quarter, Group, TVC_t, VUL_t, Sum_TVC_t, Sum_TVC_t_1, 
										  Sum_VUL_t, Sum_VUL_t_1, Sum_P_y_lm, Sum_P_y_glm, Sum_PI_start, Sum_PI_end, Sum_PI_gl_start, Sum_PI_gl_end, Sum_P_y_fix_lm )]
				  
			# Calculate ratio of revenue included into the index construction after dropping balanced groups / selection criteria in PI.pre
			PI_fix_no_nan[, ':=' (Ratio_TVC = sum(VUL_t)/Sum_TVC_t, Ratio_VUL = sum(VUL_t)/Sum_VUL_t), by = quarter]

			PI_fix_agg[, `:=`(PfixT_TVC = exp(PfixT_TVC), PfixT_VUL = exp(PfixT_VUL) )] # logged, can use sum, then unlog, for products 
			
			# Create yes indices, normalize q1 to 1, then multiply each PT_yes_TVC by the last term in the equation
			for (i in 1:nrow(PI_fix_agg)) {
				PI_fix_agg[i, `:=`(PfixT_yes_TVC = ifelse(quarter == 1, 1.0, 
							   PI_fix_agg[i, PfixT_TVC]*PI_fix_agg[i-1, PfixT_yes_TVC]),
							   PfixT_yes_VUL = ifelse(quarter == 1, 1.0, 
							   PI_fix_agg[i, PfixT_VUL]*PI_fix_agg[i-1, PfixT_yes_VUL]) )]

			}
			
			setkey(PI_fix_agg,quarter)
			setkey(PI_fix_no_nan, quarter)
		
		}
		

		# Merge original file that contained revenues with indices
		if (nrow(PI_l_no_nan)>0){
			PI_c <- unique(PI_l_no_nan)[PI_agg, nomatch = 0]
		}
			if (nrow(PI_l_no_nan)==0){
				PI_l_no_nan[,':='(PT_yes_TVC=NA,PT_yes_VUL=NULL)]
				PI_c <- PI_l_no_nan 
			}	
		if (nrow(PI_gl_no_nan)>0){
			PI_gl_c <- unique(PI_gl_no_nan)[PI_gl_agg, nomatch = 0]
		}
			if (nrow(PI_l_no_nan)==0){
				PI_gl_no_nan[,':='(PglT_yes_TVC=NA,PglT_yes_VUL=NULL)]
				PI_gl_c <- PI_gl_no_nan 
			}
		if (nrow(PI_fix_no_nan)>0) {
			PI_fix_c <- unique(PI_fix_no_nan)[PI_fix_agg, nomatch = 0]
		}
		if (nrow(PI_fix_no_nan)==0) {
			PI_fix_no_nan[,`:=`(PfixT_yes_TVC=NA,PfixT_yes_VUL=NA)]
			PI_fix_c <- PI_fix_no_nan
		}
		
		setkey(PI_c,quarter)
		setkey(PI_gl_c,quarter)
		setkey(PI_fix_c,quarter)
		
		# Right outer join to keep PI_c (Laspeyres) as master file
		# Note order is crucial to prevent empty PI_fix_c from overwriting non-empty PI_c

			PI_c <- PI_fix_c[PI_c]	
			PI_c <- PI_gl_c[PI_c]
			
			PI_c[, c("PT_TVC", "PT_VUL", "PglT_TVC","PglT_VUL","PfixT_TVC","PfixT_VUL", "TVC_t", "VUL_t"):= NULL]
			setnames(PI_c, c("Sum_TVC_t", "Sum_TVC_t_1", 
							 "Sum_VUL_t", "Sum_VUL_t_1"), 
					 c("TVC_t", "TVC_t_1", "VUL_t", "VUL_t_1"))
			rm(PI_no_nan)
		
		return(PI_c)
    }
}

# Laspeyres 
Calc.PI_L <- function(PI) {
    
    if ( !is.null(PI) ) {  
    
		# Create year variables 	
		PI[(quarter>0 & quarter<5), Year := 2006]
		PI[(quarter>4 & quarter<9), Year := 2007]
		PI[(quarter>8 & quarter<13), Year := 2008]
		PI[(quarter>12 & quarter<17), Year := 2009]
		PI[(quarter>16 & quarter<21), Year := 2010]
		PI[(quarter>20 & quarter<25), Year := 2011]
		PI[(quarter>24 & quarter<29), Year := 2012]
		PI[(quarter>28 & quarter<33), Year := 2013]	
		PI[(quarter>32 & quarter<37), Year := 2014]
		PI[(quarter>36 & quarter<41), Year := 2015]	
		
		setkey(PI, Group, Year)
		# Sum by year and Group, take mean of component indices, take lag, then use as denominator 
		New_vp <- PI[, list(V_TVC_t = sum(TVC_t),
					 V_VUL_t = sum(VUL_t), P_bar = mean(P_t, na.rm = TRUE), P_gl_bar = mean(P_gl_t, na.rm = T), P_fix_bar = mean(P_fix_t, na.rm = T)), 
					 by = key(PI)]
		# Sum by year only 
		New_vp[, `:=`(G_TVC_t = sum(V_TVC_t),
					  G_VUL_t = sum(V_VUL_t)), by = Year]
		# Weight/share of each Group 
		New_vp[, `:=`(W_TVC = V_TVC_t/G_TVC_t,
					  W_VUL = V_VUL_t/G_VUL_t)]
		# Lag of weights by Group (one year lag)
		New_vp[, `:=`(W_TVC_lag = c(0, W_TVC[-.N]),
					W_VUL_lag = c(0, W_VUL[-.N]), P_bar_lag = c(0, P_bar[-.N]), P_gl_bar_lag = c(0, P_gl_bar[-.N]), P_fix_bar_lag = c(0, P_fix_bar[-.N])),
			   by = Group]
		# For 2006, use that year's weights
		New_vp[Year == 2006, `:=` (W_TVC_lag = W_TVC,
								   W_VUL_lag = W_VUL, P_bar_lag = P_bar, P_gl_bar_lag = P_gl_bar, P_fix_bar_lag = P_fix_bar )]

		# Drop some variables)
		New_vp[,c("V_TVC_t", "G_TVC_t",
				   "V_VUL_t", "G_VUL_t"):= NULL]

		setcolorder(New_vp, c("Group", "Year", "P_bar", 
							  "P_bar_lag", "P_gl_bar", 
							  "P_gl_bar_lag", "P_fix_bar", 
							  "P_fix_bar_lag", "W_TVC", "W_TVC_lag",
							  "W_VUL", "W_VUL_lag"))
							 
		# 1. Laspeyres construction 
		PI_complete <- PI[P_t != "NA",] # Keep only if index was calculated 
		# Keep only Groups that are appear in entire sample 
		PI_complete[,max := max(quarter), by = Group]
		PI_complete <- PI_complete[max == 40, ]
		PI_complete[, max:=NULL]	
		
		setkey(PI_complete, Group, Year)
		PI_complete[, Groupkey := .GRP, by = key(PI_complete)]
		PI_complete[, num := .N, by = Groupkey]

		### In so doing, we can accomodate selection condition
	    # Selects only Groups present for all 4 quarters within a year 
		PI_complete <- PI_complete[(num == 4),]
		PI_complete[, c("Groupkey", "num"):= NULL]

		PI_merged <- PI_complete[New_vp, nomatch = 0]
		rm(PI_complete)

		# Group weights lagged one year and fixed within year, note it's a weighted arithmetic average, not geometric
		PI_merged[, `:=`(P_L_TVC = W_TVC_lag*(P_t/P_bar_lag), 
						 P_L_VUL = W_VUL_lag*(P_t/P_bar_lag) )]

		PI_L <- PI_merged[, list(P_L_TVC = sum(P_L_TVC), 
								 P_L_VUL = sum(P_L_VUL) ), by = quarter]
		setkey(PI_L,quarter)

		PI_end <- PI_merged[quarter %in% c(8,12,16,20,24,28,32,36,40),
				list(quarter, Group, P_t, P_bar, W_VUL, W_TVC)]
		rm(PI_merged)
		
		# Take out end values for each year 
		PI_end[, `:=`(P_L_TVC = W_TVC*(P_t/P_bar), 
					  P_L_VUL = W_VUL*(P_t/P_bar) )]
		
		PI_e <- PI_end[, list(P_e_TVC =sum(P_L_TVC), 
							  P_e_VUL = sum(P_L_VUL) ), by = quarter]
		setkey(PI_e,quarter)
		rm(PI_end)
		
		# 2. Geometric Laspeyres construction 
		PI_gl_complete <- PI[P_gl_t != "NA",] # Keep only if index was calculated 	
		# Keep only Groups that are appear in entire sample 
		PI_gl_complete[,max := max(quarter), by = Group]
		PI_gl_complete <- PI_gl_complete[max == 40, ]
		PI_gl_complete[, max:=NULL]	

		setkey(PI_gl_complete, Group, Year)
		PI_gl_complete[, Groupkey := .GRP, by = key(PI_gl_complete)]
		PI_gl_complete[, num := .N, by = Groupkey]	

		PI_gl_complete <- PI_gl_complete[(num == 4),]
		PI_gl_complete[, c("Groupkey", "num"):= NULL]
		
		PI_gl_merged <- PI_gl_complete[New_vp, nomatch = 0]
		rm(PI_gl_complete)

		PI_gl_merged[, `:=`(P_gl_TVC = W_TVC_lag*(P_gl_t/P_gl_bar_lag),
						 P_gl_VUL = W_VUL_lag*(P_gl_t/P_gl_bar_lag) )]

		PI_gl_L <- PI_gl_merged[, list(P_gl_TVC = sum(P_gl_TVC),
								 P_gl_VUL = sum(P_gl_VUL) ), by = quarter]
		setkey(PI_gl_L,quarter)

		PI_gl_end <- PI_gl_merged[quarter %in% c(8,12,16,20,24,28,32,36,40),
				list(quarter, Group, P_gl_t, P_gl_bar, W_VUL, W_TVC)]
		rm(PI_gl_merged)
		
		# Take out end values for each year 
		PI_gl_end[, `:=`(P_gl_TVC = W_TVC*(P_gl_t/P_gl_bar),
					  P_gl_VUL = W_VUL*(P_gl_t/P_gl_bar) )]
		
		PI_gl_e <- PI_gl_end[, list(P_gl_e_TVC = sum(P_gl_TVC),
							  P_gl_e_VUL = sum(P_gl_VUL) ), by = quarter]
		setkey(PI_gl_e,quarter)
		rm(PI_gl_end)
		
		# 3. Fixed index construction
		
		PI_fix_complete <- PI[P_fix_t != "NA",] # Keep only if index was calculated 	
		# Keep only Groups that are appear in entire sample  
		PI_fix_complete[,max := max(quarter), by = Group]
		PI_fix_complete <- PI_fix_complete[max == 40, ]
		PI_fix_complete[, max:=NULL]	

		setkey(PI_fix_complete, Group, Year)
		PI_fix_complete[, Groupkey := .GRP, by = key(PI_fix_complete)]
		PI_fix_complete[, num := .N, by = Groupkey]	

		PI_fix_complete <- PI_fix_complete[(num == 4),]
		PI_fix_complete[, c("Groupkey", "num"):= NULL]	
		
		PI_fix_merged <- PI_fix_complete[New_vp, nomatch = 0]
		rm(PI_fix_complete)		
		
		PI_fix_merged[, `:=`(P_fix_TVC = W_TVC_lag*(P_fix_t/P_fix_bar_lag),
						 P_fix_VUL = W_VUL_lag*(P_fix_t/P_fix_bar_lag) )]

		PI_fix_L <- PI_fix_merged[, list(P_fix_TVC = sum(P_fix_TVC),
								 P_fix_VUL = sum(P_fix_VUL) ), by = quarter]
		setkey(PI_fix_L,quarter)

		PI_fix_end <- PI_fix_merged[quarter %in% c(8,12,16,20,24,28,32,36,40),
				list(quarter, Group, P_fix_t, P_fix_bar, W_VUL, W_TVC)]
		rm(PI_fix_merged)
		
		# Take out end values for each year 
		PI_fix_end[, `:=`(P_fix_TVC = W_TVC*(P_fix_t/P_fix_bar),
					  P_fix_VUL = W_VUL*(P_fix_t/P_fix_bar) )]
		
		PI_fix_e <- PI_fix_end[, list(P_fix_e_TVC = sum(P_fix_TVC), 
							  P_fix_e_VUL = sum(P_fix_VUL) ), by = quarter]
		setkey(PI_fix_e,quarter)
		rm(PI_fix_end)	
		
		# Merge end values with current values 
		PI_L <- PI_e[PI_L, nomatch = NA]
		rm(PI_e)
		
		# Merge other construction methods
		# Right outer join to keep PI_L (Laspeyres) as master file
		PI_L <- PI_gl_L[PI_L]
		PI_L <- PI_gl_e[PI_L]
		PI_L <- PI_fix_L[PI_L]
		PI_L <- PI_fix_e[PI_L]	

		return(PI_L)
	
    }
}

# Take values and chain them all together to calculate indices 
	# Yes indices already chained, just working on the nsub indices 
Calc.Final <- function(PI_L, PI_Torn) {

    if ( nrow(PI_L)>0 & nrow(PI_Torn)>0 ) {
		for (i in 1:8) {

			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 1, 1.0, 
			PI_L[i-1, PL_TVC_nsub]*(PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC])),

						 PL_VUL_nsub = ifelse(quarter == 1, 1.0,
			PI_L[i-1, PL_VUL_nsub]*(PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL])), 
			
						PGL_TVC_nsub = ifelse(quarter == 1, 1.0, 
			PI_L[i-1, PGL_TVC_nsub]*(PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC])), 
			
						PGL_VUL_nsub = ifelse(quarter == 1, 1.0, 
			PI_L[i-1, PGL_VUL_nsub]*(PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL])),		
			
						PFL_TVC_nsub = ifelse(quarter == 1, 1.0, 
			PI_L[i-1, PFL_TVC_nsub]*(PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC])), 
			
						PFL_VUL_nsub = ifelse(quarter == 1, 1.0, 
			PI_L[i-1, PFL_VUL_nsub]*(PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL])) )]
				
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 1, 1.0, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]), 

							PFL_nsub = ifelse(quarter == 1, 1.0, 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm]) )]	

		}

		for (i in 9:12) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 9, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 8, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 9, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 8, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 9, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 8, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 9, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 8, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 9, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 8, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 9, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 8, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 9, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),

							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]	
		
		
		
		}


		for (i in 13:16) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 13, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 12, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 13, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 12, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 13, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 12, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 13, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 12, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 13, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 12, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 13, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 12, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 13, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),  
			
							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]
			
		}


		for (i in 17:20) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 17, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 16, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 17, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 16, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 17, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 16, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 17, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 16, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 17, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 16, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 17, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 16, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 17, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),

							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]			
		}

		for (i in 21:24) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 21, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 20, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 21, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 20, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 21, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 20, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 21, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 20, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 21, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 20, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 21, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 20, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 21, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),

							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]		
		}

		for (i in 25:28) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 25, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 24, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 25, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 24, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 25, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 24, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 25, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 24, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 25, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 24, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 25, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 24, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 25, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),

							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]	
		}

		for (i in 29:32) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 29, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 28, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 29, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 28, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 29, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 28, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 29, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 28, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 29, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 28, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 29, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 28, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 29, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),

							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]	
		}
		
		for (i in 33:36) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 33, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 32, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 33, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 32, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 33, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 32, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 33, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 32, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 33, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 32, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 33, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 32, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 33, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),

							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]	
		}

		for (i in 37:40) {
			PI_L[i, `:=`(PL_TVC_nsub = ifelse(quarter == 37, 
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[quarter == 36, P_e_TVC],
			PI_L[i-1, PL_TVC_nsub]*PI_L[i, P_L_TVC]/PI_L[i-1, P_L_TVC]),


				PL_VUL_nsub = ifelse(quarter == 37, 
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[quarter == 36, P_e_VUL],
			PI_L[i-1, PL_VUL_nsub]*PI_L[i, P_L_VUL]/PI_L[i-1, P_L_VUL]), 
			
				PGL_TVC_nsub = ifelse(quarter == 37, 
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[quarter == 36, P_gl_e_TVC],
			PI_L[i-1, PGL_TVC_nsub]*PI_L[i, P_gl_TVC]/PI_L[i-1, P_gl_TVC]),


				PGL_VUL_nsub = ifelse(quarter == 37, 
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[quarter == 36, P_gl_e_VUL],
			PI_L[i-1, PGL_VUL_nsub]*PI_L[i, P_gl_VUL]/PI_L[i-1, P_gl_VUL]), 

				PFL_TVC_nsub = ifelse(quarter == 37, 
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[quarter == 36, P_fix_e_TVC],
			PI_L[i-1, PFL_TVC_nsub]*PI_L[i, P_fix_TVC]/PI_L[i-1, P_fix_TVC]),


				PFL_VUL_nsub = ifelse(quarter == 37, 
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[quarter == 36, P_fix_e_VUL],
			PI_L[i-1, PFL_VUL_nsub]*PI_L[i, P_fix_VUL]/PI_L[i-1, P_fix_VUL]) )]
			
			PI_Torn[i, `:=` (PL_nsub = ifelse(quarter == 37, 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_PI_start]/PI_Torn[i, Sum_PI_end], 
			PI_Torn[i-1, PL_nsub]*PI_Torn[i, Sum_P_y_lm]/PI_Torn[i-1, Sum_P_y_lm]),

							PFL_nsub = 
			PI_Torn[i-1, PFL_nsub]*PI_Torn[i, Sum_P_y_fix_lm]/PI_Torn[i-1, Sum_P_y_fix_lm] )]	
		}	

		PI_L[,c("P_L_TVC", "P_L_VUL", 
				"P_e_TVC", "P_e_VUL"):= NULL]
		
		# Righter outer join using PI_L as master file
		PI_table <- PI_Torn[PI_L]
		# PI_table <- PI_Torn[PI_L, nomatch = 0]
		rm(PI_Torn, PI_L)
		
		PI_table <- PI_table[,.(quarter,Group,TVC_t,TVC_t_1,VUL_t,VUL_t_1,Ratio_TVC,Ratio_VUL,
					PT_yes_TVC,PT_yes_VUL,PglT_yes_TVC,PglT_yes_VUL,PfixT_yes_TVC,PfixT_yes_VUL,
					PL_nsub,PFL_nsub,PL_TVC_nsub,PL_VUL_nsub,PGL_TVC_nsub,PGL_VUL_nsub,PFL_TVC_nsub,PFL_VUL_nsub)]
		
		return(PI_table)
    }
	if ( nrow(PI_L)==0 | nrow(PI_Torn)==0 ) {
		PI_table <- data.table()
		return(PI_table)
	}
}














